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ABSTRACT 

In this paper we attempt to investigate the nature of the first gravitational wave (GW) signal 
to be detected by pulsar timing arrays (PTAs); will it be an individual, resolved supermassive 
black hole binary (SBHB), or a stochastic background made by the superposition of GWs 
produced by an ensemble of SBHBs? To address this issue, we analyse a broad set of simula¬ 
tions of the cosmological population of SBHBs, that cover the entire parameter space allowed 
by current electromagnetic observations in an unbiased way. For each simulation, we con¬ 
struct the expected GW signal and identify the loudest individual sources. We then employ 
appropriate detection statistics to evaluate the relative probability of detecting each type of 
source as a function of time for a variety of PTAs; we consider the current International PTA, 
and speculate into the era of the Square Kilometre Array. The main properties of the first 
detectable individual SBHBs are also investigated. Contrary to previous work, we cast our 
results in terms of the detection probability (DP), since the commonly adopted criterion based 
on a signal-to-noise ratio threshold is statistic-dependent and may result in misleading conclu¬ 
sions for the statistics adopted here. Our results confirm quantitatively that a stochastic signal 
is more likely to be detected first (with between 75 to 93 per cent probability, depending on 
the array), but the DP of single-sources is not negligible. Our framework is very flexible and 
can be easily extended to more realistic arrays and to signal models including environmental 
coupling and SBHB eccentricity. 
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1 INTRODUCTION 

Our current knowledge of the Universe is almost entirely based 
on electromagnetic observations (from radio waves, up to gamma 
rays). A direct gravitational wave (GW) detection will open 
an alternative and complementary perspective on the Cosmos 
( [Sathyaprakash & Schutz|2009^ . According to General Relativity, 
accelerating masses with a mass quadrupole moment varying in 
time generate GWs, and our Universe has plenty of systems fulfill¬ 
ing such a requirement ( |Wald|198^|Maggiore|200^ ; from compact 
binaries (emitting GWs while inspiraling towards one another un¬ 
til coalescence, Peters & Mathews|1963|>, to asymmetric spinning 


compact stars (|Van Den Broeck|2005 i, and catastrophic supernovae 


explosions jYakunin et al.|2010) , just to mention a few. 

A direct GW detection is a long-standing challenge in experi¬ 
mental physics, dating back to the ’60s. To date, two different ap¬ 
proaches to the problem have polarised the attention of the scien- 
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tific community, bringing the promise of the first direct GW de¬ 
tection. One technique is ground based laser interferometry; the 
second generation of ground based detectors (GBDs), including 
advanced LIGO (LIGO Scientific Collaborationj 2009[ l, advanced 
Virgo ( [Accadia et al.|201^ and KAGRA ( |Somiya|2012| l, will soon 
be online, and is expected to observe tens of coalescing neutron 
star or stellar-mass black hole binaries ( [Abadie et al.||2010t . The 
other technique consists of timing an ensemble of ultra stable mil¬ 
lisecond pulsars, forming what is commonly referred to as a pulsar 
timing array (PTA). The times of arrival (ToAs) of radio pulses 
emitted by galactic milli-second pulsars are collected by 100-m 
class radio telescopes around the globe. Since these pulses are so 
extremely regular, a GW passing by would introduce irregularities 
in their ToAs, a detectable fingerprint that can be measured, pro¬ 
vided 100 ns timing precision is achieved on a large number of 
pulsars ( |Sazhin|[T978[|Detweiler|197^|Hellings & Downs|1983^ . 
Regardless of which approach, either GBDs or PTAs, achieves a 
GW detection first, the two independent detections are necessary, 
because they target orthogonal classes of sources, therefore provid- 
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ing complementary information about the Universe. The invaluable 
scientific promises of GW astronomy are potentially so revolution¬ 
ary that the design of a third generation of GBDs -the Einstein 
telescope (ET)- is now being investigated ( [Punturo et al.||20T0| ), 
and the European Space Agency (ESA) has selected The Gravita¬ 
tional Universe science theme ( |eLISA Consortium]|2013| l for the 
L3 launch slot, with eLISA -the evolved laser interferometer space 
antenna- put forward as strawman mission design, now scheduled 
for 2034. 

Within the LIGO Scientific CollaboratiorQ different groups 
have formed in order to specialise in the detection of the four dif¬ 
ferent types of expected GW signals: (i) a stochastic GW back¬ 


ground (GWB,|LIGO Scientific Collaboration & Virgo Collabora- 
tion|2014bt, (ii) continuous waves (|LIGO Scientific Collaboration 


& Virgo Collaboration|2014at, (iii) bursts jLIGO Scientific Collab¬ 


oration & Virgo Collaboration|2014cj), and (iv) compact binary co- 


alescences (aka CBC, LIGO Scientific Collaboration & Virgo Col-| 
|laboration|2013^ . A GWB is the superposition of numerous similar 
GW sources, adding up in an incoherent fashion, and can have ei¬ 
ther an astrophysical origin, e.g. coming from a large population of 
a specific class of astrophysical objects ( |Rosado|201 l[[^gimbau| 
|2011| >, or cosmological, e.g. coming from physical processes such 
as inflation or phase transitions in the early Universe (|Allen|1996| 
|Maggior^|2000^ . Continuous waves are quasi-sinusoidal signals, 
expected to be produced by rotating neutron stars. The burst search 
targets short GW transients (< 1 second) with limited assumptions 
on the waveform. On the contrary, the CBC search relies on accu¬ 
rate models for the waveform emitted by coalescing binaries, and 
seeks for patterns in the detector’s data that match such models. 

PTAs are following a similar path, but maintain a flexible 
structure. There are three independent PTA consortia: the Euro¬ 
pean PTA (EPTA, [Eerdman et al.| |2010t , the Parkes PTA (PPTA, 
[Manchester et al.|2013| >, and NANOGrav ( jjenet et al.|2009| >. The 
three groups join forces and combine data under the aegis of the 
International PTA (IPTA, [Hobbs et al.|2010| l, which is formally a 
consortium of consortia. The main target of PTA campaigns are the 
GWs emitted by supermassive black hole (SBH) binaries (SBHBs; 
see, e.g., [Sesana et al.|2008^ , although processes in the early Uni¬ 
verse may also produce GWs in the same frequency band, like cos¬ 
mic strings ( Sanidas^et^ay |^12} or the cosmological QCD phase 
transition ( [Caprini et al.|2010{ ! 

On the basis of the hierarchical paradigm of structure forma¬ 
tion -according to which galaxies we see today underwent a series 
of merger events ( [White & Rees|1978^ - and on the observational 
fact that SBHs are ubiquitous in galaxy centres (e.g., [Magorrianj 
[et al.|19 98t, following galaxy mergers, a vast number of SBHBs is 
expected to form along the cosmic history ( [Begelman et al.|1980[ (. 
Although there is plenty of observational evidence of the existence 
of SBH pairs separated by hundreds-to-thousands of parsecs, de¬ 
tecting parsec and sub-parsec scale SBHBs driven by GWs has 
been proven challenging, and only candidate systems supported 
by circumstantial evidence have been proposed to date (see [Dotti| 
[et al.[[20T2l for an extensive review on the topic). The recent ac¬ 
cess to vast time domain optical data set resulted in the discovery 
of a number of quasars showing periodic variability ([Graham et al.| 
[2015[ [Liu et aL][2015[ , but the SBHB nature of these objects is 
difficult to prove. In any case, a large population of SBHBs is an 
inevitable prediction of current galaxy formation models, and their 
GW signals must have travelled to our galaxy imprinting their sig- 
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nature in the pulsars’ ToAs. The incoherent superposition of many 
signals will possibly result in a stochastic GWB jRajagopal & R o-[ 
[mani|19^ , but particularly massive and/or nearby SBHBs might 
dominate the signal budget and show up as resolvable continuous 
waves jSesana et al.|2009 1 in the data. Moreover, bursts from partic¬ 


ularly eccentric binaries ( Amaro-Seoane et al.|2010[ Einn & Lom- 
men|2010 1 or from the GW memory effect ^Eavata|2009 jPshirkov 


et al.|2010l[van Haasteren & Levin|2010^ are also possible. Conse¬ 


quently, several IPTA projects have been established, targeting dif¬ 
ferent types of signals. Unlike GBDs, there are no CBC searches 
within the PTA, since SBHBs are observed long before coales¬ 
cence, during their slow, adiabatic inspiral. For PTAs, the first de¬ 
tection will probably be either a continuous wave signal from a 
SBHB (we will refer to such signal as a ’single source’ from now 
on) or a GWB, and we focus on these two classes of sources in this 
study. 

A GWB is commonly characterised by an amplitude and a 
spectral slope ( [Jenet et^[2005[ . Detecting a GWB, with ampli¬ 
tude and slope consistent with the expected signal generated by an 
ensemble of SBHBs, would prove the validity of the models that 
predict the existence of tight SBHBs, and, in a broader context, 
would support our current understanding of the formation and evo¬ 
lution of galaxies. Moreover, the detection of a GWB could put 
constraints on alternative theories of gravity, or prove their validity 
against General Relativity ( [Chamberlin & Siemens[|2012^ . How¬ 
ever, an ensemble of SBHBs could also produce a spectrum differ¬ 
ent than a simple power law (Sesana et al.|2004l Kocsis & Sesana[ 


|20TT||Sesma|2013[[McWilliams et al.|2014[[Ravi et al.|2014t , and 

other GW sources could also lead to different measurable GWBs; 
detecting the actual shape of the GWB would help discriminate be¬ 
tween these different models ( [Sampson et al.|20T5) l. 

The detection of a single SBHB ([Yardley et al.|2010[[Corbin &[ 
[Comish|2010|[Ellis|2013[|Wang et al.|2014) would allow for some 

further appealing prospects. First, by detecting a purely stochastic 
GWB one cannot, in principle, be sure of the kind of sources that 
originated it; however, detecting a single SBHB would be a more 
direct proof of the existence of tight binaries. Detecting a single 
source would allow us to measure the characteristics of the binary, 
for example its luminosity distance, mass, orbital period, and sky 
location, opening fascinating prospects for multimessenger astron 


omy (Sesana et al. [2012)[Tanaka et al.|20T2{[Burke-Spolaor|2013[ 


[Rosado & Sesana|2014[ l. In some cases, it could also help to im¬ 

prove our measurement of the distance to the pulsars in the array 
( [Lee et al.|201 l[ l. 

Although it is commonly believed that a GWB detection will 
likely precede the identification of any single source, this statement 
has never been properly quantified in the literature. The aim of this 
paper is precisely to answer the question: what type of signal will 
more likely be detected first by a PTA, and with what probability? 
To tackle the problem, we analyse a large set of simulated reali¬ 
sations of the ensemble of SBHBs, and compare the detectability 
of the two types of signals as a function of time in each of them. 
The realisations are produced following the prescriptions of [Sesan^ 
( [2013[ , employing, for this exploratory study, the simplifying as¬ 
sumption of circular, GW-driven binaries. 

The outline of the paper is as follows. In Section [^ we de¬ 
scribe the simulations of the SBHB cosmic population and present 
the necessary mathematical tools to evaluate the detectability of the 
GWB and of single SBHBs; finally, we describe the properties of 
the pulsar arrays assumed as GW detectors. In Sectionj^ the detec¬ 
tion probability of the GWB is compared to that of single SBHBs, 
assuming the IPTA and two configurations of the Square Kilometre 
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Array (SKA); we also study the signal-to-noise ratio (S/N) of the 
two types of signals, and comment on the implications of evaluating 
their detectability using the S/N. We then investigate the properties 
of the first detectable single binaries, including the shift between 
Earth and pulsar terms of the signal. In Section|^we comment on 
some of the main assumptions adopted and caveats of the paper, 
and summarise our findings. Additionally, in Appendix [A| we in¬ 
clude a derivation of the formulas used to evaluate the detection of 
the GWB. 


2 GW SIGNAL MODELS AND DETECTION STATISTICS 

2.1 Observationally-based GW simulations 

To construct the relevant SBHB population emitting in the PTA 
frequency band, we exploit the observationally-based approach put 
forward by |Sesan^ ( |2013 ^. We assume for simplicity that SBHBs 
are all circular and GW-driven (see |Sesana|2013| for a discussion 
of the impact of binary-environment coupling and eccentricity). 

In practice, the SBHB population depends on four ingredients: 


SBH. This mass can be accreted with a different timing with re¬ 
spect to the SBH binary merger, and in different amount on the two 
SBHs, following the scheme described in Section 2.2 of |Sesana| 
|etaL] ( |200^ . 

We combine the different ways to populate the merging galax¬ 
ies with SBHs together with the galaxy merger rates to obtain 
23400 different SBHB merger rates, consistent with current ob¬ 
servations of the evolution of the galaxy mass function and pair 
fractions at 2 < 1.3 and M > and with the empirical 

SBH-host relations published in the literature. We give equal credit 
to each model, and run 10 realisations of each, producing a total of 
234000 simulated universes, each one with a particular GW signal 
produced by the ensemble of SBHBs. This number of realisations 
is sufficient to place reasonable confidence levels for the expected 
amplitude according to current observational constraints. Our ap¬ 
proach is modular in nature, and it is straightforward to expand the 
range of models to include new estimates of all the quantities in¬ 
volved. 

The output of each realisation of the Universe is a list of SB¬ 
HBs including for each system: 


(i) The galaxy merger rate. The differential merger rate can be 
written as 

(Png _ z) T{z, M, q) dU 
dzdMdq Min 10 T{z,M,q) dz 

Here, 4>{M, z) = [dn/d\o^M]z is the galaxy mass function mea¬ 
sured at redshift z\ T{M,q,z) = \df /dq\M,z is the differential 
fraction of galaxies with mass M at redshift 2 paired with a sec¬ 
ondary galaxy having a mass ratio in the range q, q-\-5q\ t{z, M, q) 
is the typical merger time-scale for a galaxy pair with a given M 
and 5 at a given 2 ; and dtr/dz converts the proper time rate into 
redshift, and is given by standard cosmology. The reason for writ¬ 
ing Equation 0 is that (j) and T can be directly measured from ob¬ 
servations, whereas r can be inferred by detailed numerical simu¬ 
lations of galaxy mergers. We take five different galaxy stellar mass 
functions from the literature jBorch et al.|2006||Drory et al.|2009[ 
[Ilber t e t^|2010[|Muzzin et al.||2013[ |To mczak et al.||2014) and 

match them with the local mass function ( [Bell et al.|20Q3| ), to ob¬ 
tain five fiducial (j)z{M). We consider four studies of the evolution 


of the galaxy pair fraction ( 

Bundy et al.|20091|de Ravel et al.|20091 

[Lopez-Sanjuan et al.[2012| 

Xu et al.|2012b. Finally, we take merger 

time-scales r estimated by Kitzbichler & White| pOOSjl and |Lotz 


[etalldMot . 


(ii) The relation between SBHs and their hosts. We assign to 
each merging galaxy pair two SBHs with masses drawn from 13 
different SBH-galaxy relations found in the literaturej|Haring^ 

Rix||2004l IGultekin et al.|[^09| |Sani et ^|2011| [Graham et al. 

201 1| [Graham [2012| [Graham & Scott[[2013| [Beifiori et al.|201^ 

McConnell & Ma|2013[[iformendy & Ho|2013| >, spanning a broad 
range of uncertainty and including recent observations that cor¬ 
rected SBH estimates upwards. 

(iii) The efficiency of SBH coalescence following galaxy merg¬ 
ers. We simply assume that SBHBs coalesce efficiently follow¬ 
ing galaxy mergers, therefore bypassing the ‘last parsec prob¬ 
lem’ ( [Milosavljevic & Merritt|2003| l. As already stated, we do not 
dig into complications related to the SBHB-environment coupling 
( [Kocsis & Sesana][2011[ [Sesana[[2013| [McWilliams et ar]|2014 [ 
[Ravi et al.|2014^ , and assume that all systems are circular and GW- 
driven. 

(iv) When and how accretion is triggered during a merger event. 
We allow, during the merger, for some amount of accretion on each 


(i) A4, the proper chirp mass, which is related to the individual 
SBH masses mi and m 2 by A4 = [mim 2 ]®^®/[mi -|- m 2 ]^^®, 

(ii) /, the observed (redshifted) GW frequency, 

(iii) 2 , the observed redshift. 

Erom the redshift, the comoving distance to a binary is calculated 
by assuming a A cold dark matter universe, 

r=-^(\Q.rr,[l + zf + nA~^^'^dz, ( 2 ) 

iiioohj Jq 

where c is the speed of light, Hiqq = lOOkms”^ Mpc“^, and the 
cosmological parameters adopted are 

(/i, Ha) = (0.7,0.3,0.7). (3) 

To determine the strength of the GW signal detected on Earth, other 
parameters need to be specified, namely: 

• longitudinal and latitudinal spherical coordinates de¬ 
scribing the SBHB location in the sky, 

• L, the binary inclination angle with respect to the line of sight, 

• f), the GW polarisation angle, 

• $ 0 , the initial GW phase. 

When simulating the signal from single sources, the previous quan¬ 
tities will be obtained by drawing numbers from the appropri¬ 
ate uniform distributions. On the other hand, when simulating the 
GWB signal, we will simply use the sky and polarisation average 
strain ( [Finn & Thorne|2000t , 

h = A^^[a^ + P], ( 4 ) 

where 

^^_^ g5/3_^5/3[^^[l + 2]]2/3 

c^r 

is the dimensionless amplitude of the signal, and the contributions 
from the two wave polarisations are defined by 

a = 1 cos^ t, (6) 

and 

b = —2 cos L. (7) 

One can note that, using the previous formulas, the strain of the 
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Figure 1. Characteristic GW strain versus observed GW frequency. For 
each realisation of the Universe we obtain a curve hc(f), as the sum of 
the contribution from all binaries (the thin black line shows the output of 
one particular realisation). The light-grey area contains all possible values 
of he found in the realisations, whereas the dark-grey area brackets the 
5th and 95th percentiles of the he distribution. The thick black line is the 
mean he at each frequency over all realisations. A frequency bin of size 
A/ = [30yr]“'^ has been assumed. 


GWB is independent of the binaries’ exact sky location and polari¬ 
sation angle. As discussed in the next section, this is likely to have 
only a minor impact on our results. 

If we time an ensemble of millisecond pulsars for a period 
T, the overall amplitude of an incoherent superposition of GWs 
can be described in terms of a characteristic GW strain he at each 
observed frequency, which is related to the strain of the individual 
sources via 


^2 _ ^kfk 


A/ 


( 8 ) 


where k is an index running over all sources in a given frequency 
bin of width A/ = 1/T. An overview of the simulated signals 
is presented in Figure [T] in which we show he for all the realisa¬ 
tions of the SBHB cosmic population. The light grey area spans the 
range of all possible he found in the realisations, whereas the 5th 
and 95th percentiles of the he distribution are contained in the dark 
grey area. One can notice that he reaches values above 10“^® at 
almost all frequency bins. This is because, over more than 2 x 10® 
realisations of the Universe, it is likely to find some extremely mas¬ 
sive and nearby SBHB. However, such high strain values can be 
regarded as rare outliers, whereas the region between the 5th and 
the 95th percentile of the distribution is much narrower. The black 
thick line gives the mean he over all realisations at each frequency 
bin, which is consistent with previous predictions for the amplitude 
of the GWB ([Rajagopal & Romani|[T^ |Wyithe & Loeb||200^ 


ISesana et al.|2008[|Rosado|2011[ Ravi et al.|2015 1 . One example 

of an individual realisation is also added in the figure, plotted with 
a thin black line. The size of the frequency bin is A/ = [30 yr]“^, 
which corresponds to the longest observing time we will consider 
in this work. At each individual frequency, he can be either dom¬ 
inated by a single loud source, or produced by a superposition of 
several SBHBs, each contributing a sizeable share of the GW strain. 
Consequently, for detection purposes, the signal might be either de¬ 
terministic or incoherent/stochastic in nature. We now turn to the 
description of the detection strategies adopted in the two cases, 
which is necessary to assess which kind of signal will likely be 
detected first. 


2.2 Detection of a stochastic background 

Let us assume that we have a large set of realisations of the Uni¬ 
verse, all of them with similar astrophysical propertied When 
searching for a GWB, we define our detection statistic S as the 
cross-correlation between the outputs of two detectors (two pul- 
sard. After a certain observing time, each realisation of the Uni¬ 
verse has a measurement of S. We assume that the collection of 
values of S from different realisations is a stochastic process. 

In the absence of a GWB, the cross-correlation output reflects 
the properties of the noise processes involved in the measurement. 
We assume this to be a stochastic Gaussian process with probabil¬ 
ity density function (PDF) defined by a mean hq and a standard 
deviation ao, 


Po(S) = 


1 




=e 


(9) 


We further assume that the noise in all detectors is white, stationary, 
with zero mean (/rq = 0), and uncorrelated. If, on the other hand, a 
GWB is present in the data (the same GWB in all realisations), the 
detection statistic will follow a different distribution, namely 


Pi(S) = 


1 






( 10 ) 


where the mean /ii is now larger than zero, and the standard devi¬ 
ation ai is, in general, different than ao- 

Given a certain value of S measured in one realisation, we 
claim that it may contain a GWB if S' St, where St is a pre¬ 
defined detection threshold. The integral of po(S) over all values 
of S ^ St gives the false alarm probability (FAP), 


roo 

a= po{S)dS, 

J St 


( 11 ) 


which is the probability of claiming a spurious detection in the ab¬ 
sence of a GWB. Alternatively, the integral of pi (S) over all values 
of S ^ St gives the detection probability (DP), 


7 = 


poo 

/ Pi{S)dS, 
J St 


( 12 ) 


which is the probability of claiming a true detection of the GWB 
when it is indeed present. 

Introducing the complementary error function (erfc), we can 
solve the integrals of Equations HD and l |12^ , to obtain 


^erfe 

2 \V2aoJ 


and 


7 = 


St — p,i 
x/2cri 


(13) 


(14) 


Throughout the paper we fix the FAP to the value aq = 0.001. We 
can then solve for St in Equation l |13| > and replace the result into 
Equation 01 to get 


7 b = -erfc 


v/2crcerfc ^(2aq) —/ii 




(15) 


^ Throughout this section, by ‘realisations’ we do not refer to the outputs 
of the computer simulations analysed in other sections of this paper, but to 
a hypothetical set of copies of the same Universe. 

^ We assume that the optimal way to cross-correlate many detectors is to 


combine them 


in pairs (J 


lAllen & Romanol 19991. 
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This is the quantity that we will use to evaluate the detectability of 
a GWB. Later on we will explicitly give formulas for /ri, ao, and 
(Ji, but first there is an important aspect of the statistic S that needs 
to be commented on. 

The cross-correlation S is constructed as follows, 

fT/2 /•T/2 

S= dt dt' (16) 
J-T/2 J-T/2 


where T is the observation time, Si{t) and Sj{t') are the data from 
two different pulsars, and is a filter function. The latter 

must be chosen in such a way that the DP is maximised for a fixed 
value of FAP; in other words, we adopt the Neyman-Pearson cri¬ 
terion to define our statistics. Assuming that S follows a Gaussian 
distribution, one can obtain closed formulas for fj,i, ao and cri as 
a function of Q{t, t'), replace them in Equation 1 15 i, and obtain 
the filter function that maximises 7 b. However, the maximisation 
of the DP is a non-trivial task. Instead, it is customary to construct 
statistics that are optimal in the sense of maximising a proxy for 
the DP, which is the S/N. The latter can be defined in two different 
ways, /ii/ao, or /ri/ai, with different implications. Hence, there 
are two possible ways to construct the optimal statistic: 


• Maximise S/Na = / ti / cto - We will refer to the statistic con¬ 
structed in this way as the A-statistic. 

• Maximise S/Nb = ^ijai. We will refer to this as the B- 
statistic. 


The derivation of these two statistics and their properties are dis¬ 
cussed in Appendix [A| 

If the signature of the GWB in the data is small, i.e. in the 
weak signal approximation, ao and ai are almost identical, and 
the two statistics become equivalent. This assumption is usually 
adopted in the literature, especially for GBD studies ([All en & Ro-| 
|mano|19^|Rosado|2012[|Regimbau et al.|2014| (, and allows us to 
define the signal-to-noise ratio as S/N= pi/ao ~ pxjax. Under 
this assumption, the true optimal filter for the cross-correlation is 
in fact the one that maximises the S/N. Furthermore, one can set 
ao ~ ai in Equation GD and fix the DP to a particular value 70 to 
obtain 

S/N'^ Si V 2 [erfc"^(2ao) - erfc"^( 27 o)] . (17) 

If for example we set qq = 0.001 and 70 = 0.95, S/N^ « 4.74 is 
the threshold such that a larger value of S/N ensures a PAP smaller 
than 0.1% and a DP larger than 95%, which could be considered a 
confident detection. 

However, [Siemens et al.j p013[ ( already pointed out that the 
weak signal approximation is, in the long run, bound to become 
inaccurate for PTAs; and since we want to make predictions for 
future, more sensitive arrays, the approximation would certainly be 
crude. If the presence of a signal is tangible, and therefore ao and 
ai are not equal (the latter being typically an increasing function 
of time), one cannot associate pre-defined values of ao and 70 to 
a fixed S/N threshold. Instead, one has two possible ways to define 
the S/N threshold: 


• If one adopts the A-statistic, S/Na = pi/ao- The threshold 
would then be 


s/nI = V 2 


erfc ^( 2 ao) — —erfc ^( 270 ) 
ao 


(18) 


• If one instead adopts the B-statistic, S/Nb = pijai. The 
threshold would then be 


S/Nb = V 2 


— erfc ^( 2 ao) —erfc ^( 270 ) 
ai 


(19) 


In both cases, as one can see, the threshold changes in time as 
pi and ai evolve. More specifically, for fixed DP and FAP, the 
S/N threshold increases for the A-statistic and decreases for the B- 
statistic. 

To avoid confusion, we prefer to evaluate the detectability of 
the signal in terms of the DP directly, avoiding the necessity of 
defining an S/N and a certain threshold. In the following we adopt 
the B-statistic, for reasons discussed in Appendix [A| This is given 
by 

Wb = 

‘iVijShoifkiSijk 

V ^ + Skoifkimifk) + Snoifk)] + ’ 

( 20 ) 

which is a linear combination of the cross-correlations between the 
data from pulsars i and j at discrete frequency fk, 

Sijk = Si{fk)sj{fk). ( 21 ) 

Here, the ‘ * ’ and ‘" ’ denote the complex conjugate and Fourier 
transform, respectively. The detailed derivation of Equation p0[ ( 
along with the formal definition of the cross-correlation is given 
in Appendix where we also derive the three quantities that are 
necessary in order to evaluate the DP in Equation ( [15^ , namely pi, 
ao, and cri. In the following we will define the various quantities 
introduced in Equation j20[ l. 

The summation in k is over all frequency bins between /min = 
T“ ^, and /max = [ Ai] “ ^, where the size of each bin is A/ = T~^, 
and the frequency fk is assumed at the arithmetic mean of bin 
k. The parameter At is the cadence time, i.e. the typical time 
lapsed between consecutive pulsar observations. We assume a ca¬ 
dence time of 2 weeks for all pulsars and arrays. Eor such a ca¬ 
dence, one has /max ~ 2 x 10“® Hz, but it is extremely unlikely 
to find a SBHB emitting at such a high GW frequency, and we use 
/max « 3 X 10“^ Hz instead in our computation. The observation 
time T is the duration of a PTA campaign; for example, the cur¬ 
rent IPTA has been recording To As of radio pulses for ~10yr. The 
other summation in the Equation l [20[ l accounts for all possible pul¬ 
sar pairs, 

M M 

E = EE- (22) 

z — l j>% 

where M is the number of pulsars in the array, which is, for the 
current IPTA, M = 49; when considering SKA-type arrays we 
will assume different values for M. Pi is the noise power spectrum 
of the z-th pulsar; we assume that the pulses have a certain degree 
of irregularity, which is a random Gaussian process described by a 
root mean square (rms) value af so that Pi is simply 

Pi = 2aiAt. (23) 


Typical values of a are between 500 ns and 5 ps (which correspond 
to the best and worst case within the IPTA, respectively). Tij is 
the overlap reduction function ( [Finn et al.[2009[|Thrane & Romano| 
[2013[ l, that for a PTA has the form ( Hellings & Downs[1983[ 


"ij = ipij In (7u) - + 2 + 2 '^7'. 


(24) 


where 7 ij = [1 — cos( 6 ij)]/ 2 , and dij is the relative angle between 
pulsars i and j. The term multiplying the Kronecker Delta 5ij is ir¬ 
relevant in the calculations, and is introduced just to normalise the 
overlap reduction function in such a way that F^i = 1. A full math¬ 
ematical derivation of Tij can be found in[Anholm et al.[([2009[(. 
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Sh is the one-sided power spectral density of the GW signal in the 
timing residuals, 


127r2/3’ 


(25) 


where he is the characteristic GW strain given by Equation 
Finally, Sho is the expected one-sided power spectral density of the 
GW signal, which is needed to construct the optimal statistic. A 
simplistic assumption regarding the signal would be that it follows 
a power law of known slope. 


Sho = 


A'^yi 13/3 

127r2 ■' 


(26) 


where A would be a fiducial characteristic strain amplitude, e.g. 
10“^®. With this definition, the corresponding characteristic strain 
of the GWB would be he = A[f , which has the ex¬ 

pected average frequency dependence of Equation § in the limit 
of a smooth signal. However, we found that in practice there is 
little difference in the performance of the statistic between using 
the fiducial model of Equation \26) and just setting Sho = Sh for 
each realisation of Sh ■ Hence, throughout the rest of the paper we 
assume that Sho is identical to Sh, given in Equation l |25^ . This 
choice can be intricate in practice, when applied to a real detection 
pipeline (given that some prior knowledge of the shape of the yet 
undetected signal is required), but is convenient for the purposes of 
this work. In Appendix]^ we justify this choice and comment on 
its practical implications. 

By working purely with the power in each frequency bin, Sh, 
in our treatment of the stochastic background we are effectively 
smearing out each source so that the emission is isotropic over 
the sky. In reality, each signal is a point source and so the GW 
power distribution on the sky will be anisotropic. Using an isotropic 
search when the background is anisotropic is sub-optimal and will 
have a lower DP. Techniques for searching for and characterising 
anisotropic backgrounds have been developed ^Gair et al.|[20T4l >, 
which could out-perform isotropic searches in certain regimes, al¬ 
though this has not yet been investigated fully. The DP of an op¬ 
timal anisotropic background search is likely to be lower than an 
optimal isotropic search for an isotropic background of the same 
net amplitude, due to the larger number of model parameters in 
the anisotropic case. For these reasons our results may be slightly 
overestimating the detectability of the stochastic background. How¬ 
ever, the background only deviates significantly from isotropy in 
frequency bins that contain only a small number of sources, which 
are mostly the higher frequency bins. The signal to noise ratio for 
the isotropic background search tends to be dominated by the lower 
frequency bins and so it is likely that any overestimate of the DP is 
fairly small. Nonetheless, this should be investigated further in the 
future, using simulations in which the contribution from each indi¬ 
vidual binary is separately added into the residuals for each pulsar. 


2.3 Detection of a single source 

The optimal way to search for a deterministic signal whose param¬ 
eters are unknown is to adopt a matched filter ( |Schutz|1997^ . The 
waveform of a circular binary is generally described by 7 -|- 2M 
parameters. The Earth term (a single sinusoidal GW) is fully de¬ 
scribed by 7 parameters: (h, f, 8, fi, fi), t, (ho), already introduced 
in Section [TT| and each of the M pulsar terms adds an additional 
GW frequency and phase {fi, >hi) to the list of parameters. The 
construction of an adequate detection statistic depends on the func¬ 
tional form of the template, and is different for evolving and non¬ 


evolving binaries (i.e., binaries for which fi < / or /i = /> re¬ 
spectively). For the sake of simplicity, we assume binaries with 
orbital evolution time-scales shorter than the typical pulsar-Earth 
light travel time (i.e. evolving binaries). For such systems, the pul¬ 
sar and Earth terms fall at different frequencies. One can then con¬ 
struct a simple template for the coherent combination of the Earth 
terms only, discarding possible contribution to the signal from the 
pulsar terms. In this case, it has been shown that the relevant param¬ 
eter space for the construction of the signal template can be reduced 
to three dimensions only; namely the frequency / and the sky lo¬ 
cation 6, (j>. The resulting T^e-statistic is optimal in the Neyman- 
Pearson sense ( [Jaranowski & Kr61ak|2005T|Babak & Sesana|2012t 
|Ellis et al.|2012| l. In the absence of a signal, the PDF of the Te- 
statistic follows a distribution with 4 degrees of freedom, 

Po{Fe) = Tee-^% (27) 

and if the signal is present, the PDF is a non-central yfi distribution 
with 4 degrees of freedom, 

ro TT 11/2 

Pi(J-e,p) = —7i(py^)e-^'-5^ . (28) 

P 

The function 7i (a;) is the modified Bessel function of the first kind 
of order 1, and the non-centrality parameter p is exactly equal to 
the optimal signal-to-noise ratio S/Ns (whose calculation will be 
explained later in this section). 

Let us assume for a moment that we know the intrinsic param¬ 
eters of the signal we are searching for, namely /, 9, and fi. In this 
case, the FAP can be simply calculated by integrating the PDF of 
the statistic in the absence of signal as 

poo 

ai= po{Te)dTe = [1 + Te]e~^’'. (29) 

However, if those intrinsic parameters are unknown (which is the 
most general case), one has to filter the data with a number of tem¬ 
plates N that cover the relevant parameter space of possible GW 
signals. Each template is an independent trial, and since we are now 
performing the same experiment N times, the total FAP becomes 

a = 1 - [1 - (30) 


where the index i identifies the FAP in the single trial case. The 
total FAP is therefore function of the exact choice of the number 
of templates N, which will be discussed later. As in the previous 
section, we can obtain the threshold in the statistic by choosing a 
certain value of FAP, ao- Then, introducing Equation l |29| > in l |30| >, 

aQ = 1 - [1 - [I + , (31) 

which allows us to numerically obtain the threshold Te- On the 
other hand, the DP can be calculated by numerically integrating 
Equation ( |28^ , 

poo 

7i = / Pl{Te,p)dTe 

= [ dTe. (32) 

P 

This is the probability of detecting one binary (in a particular fre¬ 
quency bin). But we are not interested in any specific binary, and, 
as a matter of fact, we should consider any potentially resolvable 
SBHB in our DP calculation. So, if 7 ^ is the DP of a single source 
in a specific frequency bin, then the total probability of detecting at 
least one single source in any frequency bin is given by 

7s = l-n[l-7i]: (33) 

i 
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where the index i includes all frequency bins between /mm and /max- 
Equation j33[ > gives the DP of a single source, in analogy to Equa¬ 
tion 0 , that was the DP of a GWB; these two equations are the 
main quantities that we need to investigate in order to compare the 
detectability of a GWB and single sources. In order to solve the 
integral in Equation l |32^ , we first need to know how to calculate 
p —S/Ns, which we now explain. 

In each realisation of the ensemble of SBHBs, we select the 
strongest binary in each frequency bin. By this we mean the binary 
that produces the largest single characteristic strain within the bin. 



where the index b runs over all the binaries within the same fre¬ 
quency bin in a particular realisation. The optimal S/N for the Earth 
term signal of a SBHB, in an array of M pulsars, is given by the 
coherent addition of the signal in each individual pulsar jSesana &| 
|Vecchio|2010| ) 

1 1/2 


S/Ns = 


X^S/N? 


(35) 


where 


S/N?= ^ ^ 


Si 


X f [aF+[sin($)-sin(<I>o)]-&Ff [cos($)-cos($o)]]^df. 
Jo 

(36) 

Here there are several quantities to be defined. The amplitude A is 
the one introduced in Equation 0 . The GW phase is 


<!> = 270 ft. 

The antenna pattern functions are 

p,+ ^ 1 [771 -Pif - [h-pif 


and 


F = 


1 -f D - Pi 


[m ■ pi][n ■ pi\ 


I+ n- Pi 

where the unitary vectors introduced are 

rh = -I- [sin((/) cos(^) — sin(5) cos((/) cos(@)]i: 
— [cos(</) cos(^) -I- sin(^) sin((/) cos{9)]y 
+ [sin(^) sin( 0 )] 2 , 


(37) 


(38) 


(39) 


(40) 


n = -I- [— sin((^) sin(5) — cos(^) cosp) cos( 6 ')]a; 

-I- [cos((^) sin(5) — cos(^) sin((/)) cos(0)]y 

-I-[cos(0 sin( 6 »)] 2 , (41) 


Cl = — sin( 0 ) cos(cj})x — sin( 6 ') sln((;/i)p — cos{d)z, (42) 

and 

Pi = sm{6i) cos{(f>i)x + sm{di) sm{(f>i)y + cos{9i)z. (43) 

In these equations, cfi and 9i are the spherical coordinates of the 
f-th pulsar. Einally, the noise spectral density is 

Si = 2/S.to'i -}- Sh.rest, (44) 
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where 2lS.tcj\ is the contribution from the pulsar’s white noise, and 


Sh 


,rest — 


/ 127r2/2 


(45) 


is an additional red noise term produced by all the other SBHBs 
present at the same frequency bin. Consequently, /ic,rest is obtained 
using Equation l[^, but the summation now includes the strain of 
all signals except for the strongest one, i.e. the one selected using 
Equation ( |34^ . 

If we neglect the evolution of the GW frequency during the 
observation time, the integral in Equation l |36| l can be easily solved 
analytically. In fact, the assumption that binaries do not evolve 
considerably (i.e. they are monochromatic signals) is a very accu¬ 
rate approximation, since the time-scales of the SBHB evolution is 
much longer than realistic observation times (a discussion on this 
can be found in Section II.A. 1 of |Sesana & Vecchio|20I0| l. Assum¬ 
ing monochromatic signals and after some algebra. Equation l |36| ) 
becomes 

= s£^^ [<I>T[l + 2sin^($o)] 

-|-cos(<1?t)[— sin(cE>j.) -|- 4sin(‘I>o)] — 4sin(c|?o)] 

-I- b‘^[Fff [iDTil + 2cos^(<l?o)] + sin(cE>j.)[cos(4>T) — 4cos(<I>o)]] 
—2abFf Ff^ [24>t sin(4>o) cos(<l>o) -I-sin(<l?T)[sin(4>T) — 2sin(<&o)] 
-1-2 cos('1>t) cos(4>o) — 2 cos(<&o)]], (46) 


where 


<E>t = 270 fT. (47) 

With the previous formulas we have all the necessary mathematical 
machinery to calculate the DP in Equation ( |33^ . 

There is only one missing ingredient that we have not yet de¬ 
fined, which is the number of templates, N. A careful computation 
of the independent number of templates in the search is a cumber¬ 
some task, and it is beyond the scope of this investigation. A simi¬ 
lar calculation has been performed by A. Petiteau in the context of 
the EPTA single source data analysis (Petiteau, private communi¬ 
cation), by means of a stochastic template bank. Petiteau covered 
the {f,9,<j)) parameter space, defining independent templates with 
minimal overlap of 0.5. He computed 4276 templates considering 
a volume in the frequency band 2 - 400 nHz and the whole sky. 
Based on this calculation, we consider here N = 10“^ templates 
as a reference number. This choice is somewhat arbitrary, but we 
checked that our results have little dependence on N across two 
order of magnitudes, in the range 10® < N < 10®. 


2.4 Pulsar arrays 

We will analyse the ensemble of SBHBs in the simulated universes 
assuming different pulsar arrays: 

(i) The current IPTA. It consists of 49 pulsars with fixed sky 
positions. The rms noise is of 0.5 ps for 15 of the pulsars, 4 ps for 
25, and 5ps for 9. 

(ii) A simulated SKAl array. The exact desigij^of SKAl is still 
under debate; furthermore, even when the detailed configuration of 
antennas and bandwidths is decided, it will not be clear how many 
low-noise millisecond pulsars will be detected over time. Eor the 
purposes of this work, it is enough to assume a fixed amount of 


http://www.skatelescope.org 
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pulsars with a certain average rms noise. In particular we assume 
that the SKAl PTA will consist of 50 pulsars with noise around 
100 ns. 

(iii) A simulated SKA2 array. The characteristics of the SKA2 
configuration are even more uncertain than those of SKAl. Again, 
we adopt a simple approach and assume that the SKA2 PTA will 
be able to monitor 200 ms pulsars with noise around 50ns. 

In the SKA cases, the sky positions of the pulsars are chosen dif¬ 
ferently (and randomly) for each realisation of the Universe. We do 
this, instead of just using one particular random distribution of pul¬ 
sars over the sky, to avoid a possible bias in the results, that could be 
provoked by a particularly convenient or inconvenient distribution 
of pulsars. 

As mentioned before, we consider GW signals in the fre¬ 
quency range [/min = r“^,/max = Af“^] in frequency bins of 
width A/ = T~^, and assume that the pulsar noise is described 
by a stationary Gaussian white process defined by a certain rms 
noise. This oversimplistic approach does not take into account two 
complications: 

(i) The presence of residual uncorrelated red noise. It is likely 
that all pulsars will show intrinsic red noise behaviour to some level 
(which has in fact already been detected in some of them); however, 
this is not always necessarily the case. We therefore assume here 
no additional red noise, and defer the investigation of its effect to a 
future pape^ 

(ii) The fitting for a pulsar timing model. This is a step inherent 
to the timing process and cannot be hampered without affecting the 
final results. Most importantly, the timing model fits out a quadratic 
function g{t) = a -\- ht + cd encoding the pulsar spin-down. This 
will partially absorb the red spectrum imprinted by the GWB in 
each individual pulsar, affecting the detection capability of the ar¬ 
ray. As shown by |Moore et al.| ( [2015) , this will likely affect the 
sensitivity in the lowest couple of frequency bins (see Figure 1 in 
that paper). 

Our main goal is not to make accurate predictions about when 
a PTA will detect the first GW signal, but rather about what kind of 
signal it will be. The answer to this second question is also some¬ 
what dependent on the treatment of the timing model, but to a much 
lesser extent. We plan to implement consistently the effect of both 
timing model and additional red noise in the future. In this first pa¬ 
per we will show results assuming the whole frequency domain, 
[r“^, Af“^], removing the lowest frequency bin, [2T’“^, Af“^], 
and removing the two lowest frequency bins, Af“^]. A 

careful implementation of the effect of the timing model is expected 
to lie between the two latter cases, i.e., removing the lowest and the 
two lowest frequency bins. 


3 RESULTS 

Having described all the relevant methods, we now turn to the dis¬ 
cussion of our main results. For each of the 234000 ensembles of 
SBHBs, we compute he, using Equation progressively increas¬ 
ing the observation time T from 1 to 30 yr in steps of 1 yr for the 
IPTA, and from 0.5 to 10 yr in steps of 0.5 yr for the SKAl and 
SKA2 arrays. For each value of T we then compute the DP of the 
GWB according to the B-statistic using Equation {TsJ and the S/N 

^ This amounts to simply add an appropriate to the noise power 

spectmm of each pulsar, yielding Pj = 2(T? At + 



Time span/ yr 



Time span/ yr 

Figure 2. Detection probability (averaged over all realisations of the ensem¬ 
ble of SBHBs) versus observing time of a GWB (green lines) and a single 
source (red lines) assuming the IPTA (upper panel) SKAl (middle panel) 
and SKA2 (lower panel). Solid lines show DPs integrated over all the fre¬ 
quency range [T“^, whereas the dashed (dotted) lines show the 

result when the lowest frequency bin (two lowest bins) are removed from 
the calculation. 

from Equation ( |A19^ . We also identify the loudest sources at each 
frequency bin and compute the DP of a single source and its S/N 
according to Equations j33[ l and j35[ l, respectively. 

3.1 Detection probability as a function of time 

Although addressing when we will detect GWs with PTAs is not the 
main goal of this work, the first natural outcome of our calculation 
is the DP of each class of sources as a function of time: 7 b for a 
GWB, and 7 s for an individual source. This is presented in Figure 
l^for the three investigated arrays. In each panel, we show the result 
for a specific PTA, considering the whole signal and progressively 
subtracting the lowest frequency bins. 

It is clear from Figure|^that the subtraction of the lowest fre¬ 
quency bins has a larger impact on the DP of the GWB than on the 
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DP of single sources. This is because there are more binaries emit¬ 
ting at lower frequencies, and thus the GWB becomes stronger, and 
more likely to be detected, at the lowest frequency bins. For single 
sources, the trend is different. On the one hand, at very low frequen¬ 
cies there are more binaries, but they are hardly resolvable, i.e. the 
GWs produced by all other binaries act as a noise that masks the 
signal of an individual source. On the other hand, at high frequen¬ 
cies, although single sources are easier to resolve (since there are 
fewer binaries emitting at similar frequencies), it is less likely to 
find one emitting strong enough GWs to be detectable. There is 
thus a particular range of frequencies (not necessarily the lowest or 
the highest frequency bins) where single sources are more easily 
detectable. 

The green curves in the upper plot of Figure]^ show the pre¬ 
dicted DP of a GWB for the IPTA. The current array has been run¬ 
ning for « 10 yr, which corresponds to a DP of « 37% assuming 
pure white noise (solid green curve). However, the DP at 10 yr goes 
down to « 5% and < 1% when the lowest and the two lowest fre¬ 
quency bins are taken out of the computation. Another way to put 
this is by looking at the 50% DP: this probability is reached in ap¬ 
proximately 2 yr, 9 yr, and 16 yr from now, depending on the cut-off 
on the low frequency bins. 

The DP for individual sources with the IPTA (red curves in 
the upper plot of Figure]^ is, on the other hand, always relatively 
small, reaching values around 10 — 20% after 20 years from now. 
This should not be taken as an indication against the development 
of single source searches. Firstly, we adopted a fairly conservative 
criterion for the detection of single sources, placing a threshold of 
FAP a = 0.1%; it could be possible that the signal looked like 
a collection of relatively faint hotspots in the sky, not classified as 
single sources in our computation, but also not similar to a GWB. A 
search for multiple single sources ( |Babak & Sesana|2012^ might be 
more efficient than a search for a GWB in such situation. Secondly, 
these predictions assume the current array without any future im¬ 
provement. Adding more pulsars and, most importantly, reducing 
the rms of individual pulsars will result in a much higher chance 
of detecting a single source, which is what we expect as we en¬ 
ter the SKA era. Indeed, this is shown in the middle and lower 
panels of Figure]^ (note that the x-axis now runs to 10 year only). 
With the full SKA PTA (described in Section [2!4) the detection of a 
GWB will be almost certain within 5 years, whereas single sources 
will be detected with a probability higher than 80% after 10 years. 
The higher probability of detecting individual sources stems from 
the lower rms expected from individual pulsars in the SKA era. 
The PTA will start to ‘hit the signal’ at much higher frequencies 
(/ > 10~® Hz), where it is also likely to see individual sources and 
not only an unresolved GWB. We notice, nonetheless, that these 
numbers are quite less optimistic that those quoted in |Janssen et al.| 
( |2015^ . This is because, in that paper, both the GWB and single 
sources are said to be detectable if the produced S/N is larger than 4. 
However, it turns out that this is not sufficiently stringent for claim¬ 
ing a confident detection of a single source. One should consider 
the results of this work as an update to those claimed in| |Janssen| 
|et al.|p015| l, keeping in mind some important caveats related to 
the possible resolution of multiple sources that we will discuss in 
Section|4] 

Let us stress once again that when calculating the DP for the 
SKA, we neglect previous IPTA data. One could achieve a more ac¬ 
curate prediction by adding millisecond pulsars to the IPTA grad¬ 
ually, until reaching the SKAl configuration. This may be consid¬ 
ered for a future work, whereas in this paper we attempt to focus on 



Figure 3. S/Ns and S/Ng* versus observed GW frequency. The red area 
contains the S/Ns values of the strongest binaries in each frequency bin, 
for all of our realisations of the SBHB population; the red line gives the 
average S/Ns over all realisations. The green area contains the cumulative 
contributions to S/Nb above each paificular frequency bin, S/Ng*, defined 
in Equation |48| ; the green line gives the mean S/Ng* at each bin. Thus, the 
left-most point of the green line gives the ensemble average S/Ng. 


the detection capabilities of the current IPTA, and compare it with 
a PTA constructed with the SKA. 


3.2 Signal-to-noise ratio 

As explained in Section]^ the probability of detecting a GW signal 
is evaluated in this paper using the DP instead of the S/N, which is 
more commonly used in the related literature. We now explore the 
properties of the S/N of the two types of signals. 

This is shown in Figure[^for the IPTA; there we compare S/Nb 
with S/Ns, calculated as explained in Sections [2.2| and [23] respec- 
tively, assuming an integration time of 30yr, and a frequency bin of 
A/ = [30yr]“^. The red area contains the S/Ns of the strongest 
binary at a given frequency bin, in each of the realisations; the red 
line is the mean S/Ns over all realisations. The green area, instead, 
contains the cumulative contribution to S/Nb (whose formula will 
be derived in Appendix]^ from frequencies above a particular fre¬ 
quency bin, i.e. 

S/N^‘ = 

r . n 1/2 

/max M M ■p2 c*2 

^ A- + [P^ + PA + [1 + ry ■ 

The green line in the figure is, again, the average over all reali¬ 
sations. Note that the values of S/Ng” at each frequency bin are 
irrelevant; it is the overall values of S/Nb (which are the left-most 
points in the green area) that are actually meaningful. On the other 
hand, at each frequency bin we have independent single sources, 
some of which might be bright enough to be resolvable. There¬ 
fore all values of S/Ns are equally relevant. We point out that, at 
/ > 10“® Hz, the typical S/N of the brightest singles sources is in¬ 
deed higher than the cumulative S/N of the GWB down to the same 
frequency, whereas it flattens out at / < 10“®Hz, where there are 
more and more sources per frequency bin, and it is difficult for a 
single SBHB to dominate the signal. 

The evolution of the S/N as a function of the observation time 
is shown in Figure]^ for IPTA, SKAl and SKA2. The S/N presents 
similar features for both classes of sources, behaving as a double 
power law; however, their detailed evolution is different, and the 
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Figure 4. S/N versus observation time, for an IPTA (top), SKAl (central), 
and SKA2 (bottom). The red area covers the values between the 5th and 
95th percentiles of S/Nj of all realisations; the red line is the average over 
realisations. The green area covers values between the 5th and 95th per¬ 
centiles of S/Ng, and the green line is the average. 


two curves end-up intersecting at some point. The S/N produced 
by a GWB follows the trend already predicted by [Siemens et al.| 
( [2013^ , which is not surprising, given that the formula they use to 
calculate the S/N is almost identical to ours (numerically, the dif¬ 
ference is negligible), in Equation ( |A19^ . One thing that is worth 
noticing is that the specific value at which the S/N flattens out for 
long observing times does only depend on the number of pulsars 
in the array. In the IPTA case, where 15 pulsars give a major con¬ 
tribution to the detection statistic (having a rms noise which is 8 
times better than the others) the turnover starts already at S/N~ 1, 
whereas in the SKA2 case, where 200 pulsars equally contribute to 
the array, the turnover starts only at S/N« 10. In practice we have 
a turnover by S/N oc where M is the number of pulsars. 

Therefore, having an array formed by many pulsars with decent 
precision might be a better strategy to detect a GWB than hav¬ 
ing few ultra-precise pulsars (a point that has also been made by 
[Siemens et al.|20r^ . 


The S/N of single sources grows initially faster than quadrat- 
ically in time, because by accessing lower and lower frequencies 
(but still higher than 10“® Hz, where SBHBs are still quite sparse) 
there is a larger chance to detect a new bright individual source. 
However, after enough observation time, when frequencies lower 
than 10“® Hz start to be accessible, individual sources with high 
S/N become less likely to be found. At that stage, the S/N of an 
individual source grows as as suggested by Equation j36| . 

By inspecting the upper plot in Figure]^ one could claim that 
the detection of a single binary is currently (by « 10 yr of ob¬ 
serving time span for the IPTA) more likely than that of a back¬ 
ground, since the average S/Ns is larger than S/Nb. Moreover, if 
one assumes that a detection is achieved as soon as the S/N sur¬ 
passes a certain threshold (for example 4, as in [Janssen et al.[20T5) >, 
one would also claim that detecting a single SBHB is more likely 
than detecting a GWB for the two SKA configurations considered. 
However, both claims would be wrong; as we mentioned in Section 
1^ the detectability should be evaluated in terms of the DP, and, 
as shown in Figure]^ the probability of detecting a GWB is in all 
cases (for all PTAs and observing times considered) larger than that 
of detecting a single SBHB. 


3.3 What will be detected first? 

Having inspected the DP of each class of sources, we can now 
tackle our original problem: which kind of signal is more likely 
to be detected first by a PTA? To answer this question we first need 
to know the probability of a GWB being detected between times t 
and t + dt\ let us call this probability pB{t)dt. The probability of 
detecting a GWB after a time t would then be 

7s(f) = [ PB{t')dt', (49) 

Jo 

where 7s is the DP of a GWB, given by Equation ( [15^ . Similarly, 
the probability ps{t)dt of a single source being detected between 
times t and t -|- dt, would be related to the DP of a single source (in 
Equation ( [^ ) by 

7s(f) = [ ps{t')dt' . (50) 

Jo 

We can numerically obtain the probability density functions pB{t) 
and ps{t) from 7s (f) and 7s (f), respectively. Once those func¬ 
tions are known, we can calculate the probability that a single 
source is detected between t and t + dt given that a GWB has 
not been detected at any time before t + dt; this would be given by 
[1 — 7s(t)]ps(t)di. Finally, we can define the function 

Ps{T)=( [1 - 7fl(t)]ps(f)dt, (51) 

Jo 

which gives the cumulative conditional probability over time of a 
single source being detected after a time T, given that a GWB has 
not been previously detected. This is the quantity that we need to 
evaluate in order to answer the central question of this paper. 

We compute Ps{T) as a function of the observation time T 
for the different arrays considered in this work; the result is shown 
in Figure]^ There we see that, if the IPTA array keeps going for 
other 20 years, there is only about a 4-8% probability that the first 
claimed detection is of a single source. Thus, the first detected PTA 
signal will most likely be a GWB or an incoherent superposition 
of multiple sources. Things will change significantly (but not dra¬ 
matically) in the SKA era. An SKAl-type PTA will have a 8-13% 
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3.4 Properties of the first detectable single source 

Even though single sources will most likely be detected after a 
GWB, they will eventually pop-up in PTA data. It is therefore inter¬ 
esting to study their properties in terms of mass, redshift and GW 
frequency. We summarise this information in Figure]^ for all the 
investigated arrays and different observation times. All frequency 
bins were taken into account when constructing this figure (even 
the lowest ones). To construct these distributions, we took all the 
realisations of the SBHB population and considered, at any obser¬ 
vation time, the brightest source in each frequency bin, weighting 
their contribution by their DP. 

The properties of the first detectable single sources do not 
depend strongly on the assumed PTA. The vast majority of the 
sources are massive, nearby binaries, clustering in the redshift 
range 0.3 — 0.7, regardless of the PTA properties. In terms of chirp 
mass, the more sensitive the array, the less massive the first de¬ 
tectable binaries would be. Therefore, although the IPTA has a 
chance to detect binaries with M > lO^'^M©, SKA2 will ob¬ 
serve many more systems with A4 ~ 10® M©. The rms noise of 
the pulsars in the array also has an impact in the frequency dis¬ 
tribution of the systems (right panels). In a putative IPTA scenario, 
the higher DP is always at the lowest frequency, despite the fact that 
more sources contribute to the signal in that range, and the inherent 
probability of having a single source sticking out is smaller. This is 
simply because the individual pulsar rms in the IPTA are not good 
enough to allow efficient detection at high frequencies, and sources 
must be extremely bright to be observed there. Conversely, in both 
SKA scenarios (central and lower rows of plots), the single source 
DP peaks around 10~®Hz, even if the observation time is extended 
to 10 years and lower frequencies are accessible. This is simply 
because the array sensitivity now allows the detection of individual 
sources also at higher frequencies, making clear that the presence of 
these systems becomes intrinsically less probable at progressively 
lower frequencies. 


3.5 Pulsar term of the first single source 


Figure 5. Cumulative probability of detecting a single resolvable source 
before a GWB versus observing time, for the IPTA (top panel), SKAl (cen¬ 
tral panel) and SKA2 (bottom panel). In each panel, solid lines show DPs 
integrated over all the frequency range whereas the dashed 

(dotted) lines show the result when the lowest bin (two lowest bins) are 
removed from the calculation. 


chance to detect a single source first, a percentage that might in¬ 
crease to about 12-25% in the SKA2 phase. Note that the afore¬ 
mentioned probabilities depend on the treatment of the lowest fre¬ 
quency bin. In particular, when all low frequency bins are consid¬ 
ered, the DP of the GWB becomes much larger than that of single 
sources, and hence the probability of detecting a single source be¬ 
fore a GWB becomes lower and saturates earlier (solid curves in 
Figure]^. To be conservative, we assume that the lowest frequency 
bin has a smaller contribution in the signal build-up, because of 
the timing model fitting; then, omitting the results in which all fre¬ 
quency bins are considered, we can confidently say that the chance 
to detect a single source before a GWB lies between 7% and 25%, 
depending on the array configuration. 


The GWs arriving at our galaxy produce a deformation of the 
space-time metric that can be observed in two ways: via the Earth 
term and via pulsar terms. The Earth term affects the ToAs of the 
pulses of all pulsars at the moment when the GWs reach Earth. On 
the other hand, a pulsar term is originated when the GWs reach any 
of the pulsars; the distorted pulses emitted by that pulsar are ob¬ 
served by our telescopes after they travel all the way from the pulsar 
to Earth. Thus, there is a delay between the Earth term and any ob¬ 
servable pulsar term, of hundreds or thousands of years, depending 
on the distance to the pulsar and its relative angle with respect to the 
GW source. Consequently, it is possible to observe GWs of a sin¬ 
gle source at two different stages of the binary’s lifetime, differing 
hundreds or thousands of years. The question we address now is, 
does the GW frequency differ significantly between the Earth and 
pulsar terms? This is an interesting topic that has consequences in 
the development of data analysis algorithms targeting continuous 
GW sources. 

Let us call fi the observed GW frequency of a pulsar term, 
and the observed GW frequency of the Earth term. We assume 
that 2 is the redshift of the GW source, which can be safely con¬ 
sidered invariant during time-scales of thousands of years. Then, at 
the time when the SBHB emitted those waves, the pulsar and Earth 
terms had frequencies /e = [1 + z\fi and /e = [1 + (using a 
similar nomenclature as in|Rosado|201 1||. The interval of time that 
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Figure 6. Properties of the individual SBHBs that are most likely to be detected with the current IPTA (upper row plots), SKAl (central row) and SKA2 
(lower row). The different columns of plots show, from left to right, the probability density function of redshift, chirp mass, and GW observed frequency of 
the individual SBHBs. The black points in the graphs on the right column are located at the arithmetic mean of the frequency bins. 




Figure 7. Histogram of the frequency shift between the pulsar and Earth terms, in units of the frequency bin T~^), where T is the observing time (assumed 
10yr in all cases), for the current IPTA (left panels), SKAl (central panels) and SKA2 (right panels). The upper plots assume a typical observed time delay 
consistent with the mean estimated pulsar distance in the current IPTA (p=; 1.5 kpc). The lower plot assumes a typical time delay consistent with the maximum 
pulsar distance in the current IPTA (j^ 6 kpc). 


the SBHB needed to evolve between the emission of those waves 
is given by 

^ 2567r8/3[GA4]5/3 ^ - [/=] ' ]■ (52) 

That interval of time would be today observed as AT = AT [1 + 


z\. Introducing this in Equation l |52[ l and rearranging, we get 

Ar256^V"[GAf[l + z]]"/" + 


f = 


5c3 


(53) 


Equation ID allows us to quantify the frequency of the pulsar term 
(/^) once the chirp mass, redshift and Earth term frequency (/^) 
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The nature of the first PTA detection 13 


of a putative resolvable source (all information we get directly from 
our Monte Carlo realisations of the Universe) are known. The only 
other quantity we need to fix is the delay AT between the arrival 
of the Earth term and the pulsar term, which depends on the typical 
distances of the observed pulsars to the Earth, that we next quantify. 

The most distant pulsar in the current IPTA is at ~ 6kpc. 
The maximum delay between pulsar and Earth terms for this pul¬ 
sar would be of ~ 12kpc/c, which is twice the light travel time 
between that pulsar and Earth. This delay would be achieved if the 
SBHB was located exactly at the opposite sky location of the pul¬ 
sar. We calculate fi of all binaries in each realisation of the Uni¬ 
verse, assuming AT = 12kpc/c, and also AT = 3kpc/c, which 
would be a typical delay for average IPTA pulsars (the mean esti¬ 
mated distance to IPTA pulsars is ~ 1.5). Then we obtain the fre¬ 
quency shift between the pulsar and the Earth term, [fi — fi]/Af, 
where A/ = T~^, and T is the observing time needed to detect 
the SBHB. 

The distribution of frequency shifts of the first detectable bi¬ 
nary for all the considered PTAs is shown in Figure]^ assuming 
AT = 3 kpc/c (upper plot) and AT =12 kpc/c (lower plot). All 
histograms are weighted by the DP of the binaries, so that binaries 
with larger probability of being detected are more representative 
than those with smaller DP. By looking at the figure one can con¬ 
clude that a significant evolution between pulsar and Earth term is 
not likely; among binaries detected in the near future with the IPTA, 
a shift of at least 1 bin has a probability of 22% for average pulsars 
(at around 1.5 kpc), and to 33% for a pulsar at 6 kpc. These percent¬ 
ages increase as we move into the SKA era; for an SKA2 arrays 
the odds of having a sizeable frequency shift between the pulsar 
and the Earth terms become higher than 50%, due to the fact that 
single sources are detected at higher frequencies, where their chirp¬ 
ing becomes progressively faster. All in all, we can conclude that 
frequency-shifting and strictly monochromatic individual sources 
are roughly equally likely to be detected. 


4 DISCUSSION AND CONCLUSIONS 

The results shown in the previous section have a number of impor¬ 
tant implications for the future design of PTA observations and data 
analysis pipelines. 

(i) We showed that a GWB is more likely to be detected than any 
resolvable single source, which is consistent with other recent pre¬ 
dictions ( [Rosado & Sesana|2014[|Ravi et al.|2015t . However, this 
result should not discourage the development of single source de¬ 
tection pipelines for at least two reasons. First, there is a 5-to-25% 
probability of detecting a single source first, which is not negli¬ 
gible; in particular, if the IPTA sensitivity will prove insufficient, 
expected timing improvement in the SKA era will result in an in¬ 
crease of the odds of detecting a single source. Second, throughout 
the paper we adopted a conservative definition of GWB (isotropic, 
stochastic, Gaussian, unpolarised, and stationary); the true signal 
produced by the superposition of GWs coming from the ensemble 
of SBHBs might well be dominated by a handful of signals, there¬ 
fore significantly departing from isotropy and/or Gaussianity. The 
development of detection algorithms targeting multiple individual 
sources j Babak & Sesana|2012[|Petiteau et al.|2013[ ( as well as cer¬ 
tain types of anisotropic signals ( jGair et ^ ]20T4) might prove to 
be a ‘more optimal’ strategy than searching for a GWB with the 
aforementioned properties. 

(ii) Predictions should generally be reported in terms of DP and 
not S/N. From a frequentist perspective, what matters is the DP at 


a fixed false alarm probability; a concept which translates into dif¬ 
ferent S/N values for different types of sources and different data 
analysis techniques. In our specific case, for fixed DP and FAP, a 
much larger S/N is required for the detection of a single source than 
for the detection of a GWB, because of the large number of tem¬ 
plates needed in the search of the former. Moreover, even within 
the same class of sources, depending on the adopted statistic, fix¬ 
ing DP and FAP may result in a time-dependent S/N threshold (in 
particular when the GWB becomes no lonrer a weak signal), which 
is the case for both our A- and B-statistic^ as discussed in Section 

m 

(iii) Nonetheless, the S/N behaviour as a function of time en¬ 
codes interesting information about source detectability. Both the 
S/N of a GWB and of a single source increase rapidly with time 
in the weak signal limit, but after the weak signal regime is aban¬ 
doned the increase becomes slower. The transition between the two 
regimes appears to happen at an S/N oc , as shown in Figure 

Therefore, adding to a PTA many pulsars with decent precision 
might be a better strategy to detect any type of signal than hav¬ 
ing few ultra-precise pulsars (a point that has also been made by 
[Siemens et al.|20r3] in the GWB case), and this issue is certainly 
worthy of further investigation. 

(iv) The first single sources to be detected will likely be ex¬ 
tremely massive {M > 10®Mq), and emitting at a GW frequency 
of «10nHz. Incidentally, for typical millisecond pulsar distances 
of ~kpc, the transition between evolving/non-evolving sources oc¬ 
curs at about that same frequency. Therefore, it is impossible to 
discard a priori waveform models and/or detection techniques. For 
example, from a frequentist perspective, both Tp and Te statistics 
are equally important tools, and no option should be discarded. 


Those results were obtained under a number of simplifying 
assumption that we now discuss. 

• We assumed circular, GW-driven binaries. Although this as¬ 
sumption might be reasonable, SBHBs observable with PTAs, are 
at the centi-to-milliparsec separations, where the coupling to the 
environment can still play a role. This fact has been pointed out 
by a number of authors jKocsis & Sesana|[2011[ [Sesana[[2013[ 


[McWilliams et al.|2014[[Ravi et al. 20141, and might be relevant for 
the results reported above. Efficient environmental coupling drives 
SBHBs faster towards coalescence; the net effect is that there are 
less systems contributing at each frequency bin, and the odds that 
a single source will dominate the signal are likely to increase. We 
therefore speculate that environmental coupling, although damag¬ 
ing in terms of GWB detection, might promote the detection of sin¬ 
gle sources. On the other hand, environmental coupling might also 
result in very eccentric binaries ( [Sesana[f2010[ [Roedig & Sesana[ 
[2012[ (. In this case, each individual system will emit in a whole 
range of GW harmonics, resulting in a more complicated signal. 
The development of data analysis pipelines for eccentric SBHBs 
has recently begun ( [Zhu et al.|2015[ l, but the DP of such systems 
requires further investigation. Regardless of the precise mechanism 


® We note that this does not apply to the statistic employed by |Siemens| 
[ef^pOB) , which is similar to our B-statistic, but with a subtle differ¬ 
ence in the definition of the null hypothesis: they assume that all pulsars 
show a common but uncorrelated red noise (due to some unknown physical 
process) whose spectral shape mimics exactly that of the GWB. Under this 
conservative assumption, <jq = and the S/N threshold corresponding to 
a given DR at a fixed FAP is indeed constant in time. A more detailed study 
of the statistic used in|Siemens et al.H2013^ can be found in|Chamberlin| 
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that could affect the actual shape of the GWB, their overall effect 
could be a lower amplitude of the signal around the lowest fre¬ 
quency bins. Thus, in practice, subtracting the contribution from 
the lowest frequency bins (which we did in previous sections in or¬ 
der to simulate the effect of the timing model) could be a sensible 
way to mimic the reduction in the signal. In such a case, the DP, 
plotted in Figure]^ would be well modelled by the dashed and dot¬ 
ted curves, that disregard the signal from the lowest and two lowest 
bins, respectively. 

• When evaluating the DP of single sources, we considered only 
the loudest binary at each frequency bin, and treated the rest of the 
signal as a stochastic GWB. However, such approximation might 
be too crude, since the rest of the signal might be dominated by just 
a few SBHBs. [Babak & Sesana| ( [20T^ showed that an array of M 
equally good pulsars can in principle resolve M/3 sources at each 
frequency. Although the IPTA data set is dominated by a few, very 
good pulsars, and the resolution of multiple sources seems unlikely, 
things will be different in the SKA era. Having a large number of 
extremely stable pulsars might realistically favour the resolution of 
multiple SBHBs at each frequency, and, in this respect, the 80% DP 
probability for singles sources with an SKA2-type array should be 
regarded as a robust lower limit. 

• When calculating the detectability of a single source, the sky 
location of the SBHBs is taken into account; one can see this in 
Equation ( |46^ , that depends on the sky coordinates of the binary 
and the pulsar. However, the DP of a GWB is insensitive to the 
particular distribution of the binaries over the sky; the function Sh, 
that encompasses the information about the signal in the DP, de¬ 
pends on the characteristic amplitude of the GWB, made up as the 
superposition of the GW strains from all SBHBs in the ensemble 
(summed in quadrature). This implies that the particular location 
of the binaries on the celestial sphere is not relevant (only the dis¬ 
tance to the binaries is), and the contribution of the ensemble to 
the cross-correlation is equivalent to that of an isotropic GWB. In 
other words, the signal of each binary is isotropically smeared out 
over the sky. This can favour the predicted detectability of a GWB 
over single source^ since the statistic employed for the detection 
of the background is optimal for an isotropic signal. If the ensem¬ 
ble of SBHBs produced a rather anisotropic signal, it could still be 
detected with a cross-correlation analysis, but the result would be 
suboptimal ( [Cornish & Sesana|2013|). In su ch a case a search for an 
anisotropic GWB i Mingarelli et al.|2013^ or a search for multiple 
sources ( jBabak & Sesana|2012 1 would be more efficient, although 
the performance of an isotropic search or an optimal anisotropic 
search when the background is anisotropic is likely to be a bit lower 
than that of an optimal isotropic search for an isotropic background, 
which could delay detection. A more sophisticated way to evaluate 
the detectability of the GWB would require the DP to be sensitive 
to the sky distribution of binaries, by considering the contribution 
of each binary to the residuals of each pulsar in the array. Given 
the large amount of simulations analysed in this work, such a task 
would presumably become computationally involved. We point out, 
nonetheless, that the DP of the GWB has most of its contribu¬ 
tion from the lowest frequency bins, where the background in most 
models can safely be considered isotropic. Significant anisotropies, 
that arise when the background is dominated by fewer binaries, are 
usually important at rather high frequencies, that do not contribute 
as much to the DP. Therefore, we expect that this possible over¬ 
estimate of the GWB DP should not be significant, and evaluating 


^ We thank Neil Cornish for pointing out this possible caveat. 


the detectability of the GWB with a more sophisticated approach 
would not change the conclusions of this work. 

• We considered ideal millisecond pulsars described by white 
noise only, ignoring the impact of fitting for a timing model. This 
issue was already commented on in previous sections, as well as the 
possibility of including unmodelled red noise in our estimates. We 
defer a more comprehensive treatment of the pulsar noise model to 
a future work. 

• We considered very idealised PTAs, assuming that the IPTA 
will simply continue as it is now, without adding any further pulsar, 
or improving their timing precision. Moreover, we kept IPTA and 
SKA separated; a realistic computation of PTA detection probabil¬ 
ities in the SKA era should take into account all previous, valuable 
IPTA data. 

Despite all these limitations and caveats, the calculation we 
presented here is the first quantitative attempt at assessing the ‘sin¬ 
gle versus background’ issue. An unresolved GWB is more likely 
to be detected, but there is a sizeable chance to see a single source 
first, and the development of the relevant detection pipelines should 
not be stopped. An extension of this work to more realistic PTAs 
and source populations including environmental coupling and ec¬ 
centricity will be crucial to direct the development of data analysis 
pipelines and to possibly bring the first GW detection with PTAs 
closer in time. 


ACKNOWLEDGEMENTS 

The authors thank the useful help and comments from Bruce Allen, 
Carsten Aulbert, Stas Babak, Ewan Barr, Neil Cornish, Justin El¬ 
lis, Ilya Mandel, Giulio Mazzolo, Antoine Petiteau, Joe Romano, 
Stephen Taylor, Yan Wang, Karl Wette, and Xingjiang Zhu. AS and 
JG are supported by the Royal Society through the University Re¬ 
search Fellow scheme. 


APPENDIX A: OPTIMAL STATISTICS EOR 
GRAVITATIONAL WAVE BACKGROUND DETECTION 

AI Cross-correlation 


Here we derive and discuss the properties of the two statistics, A 
and B (introduced in Section [2^ , that can be used to characterise 
the detection of a GWB. For a filter function of the form Q{t — 
t'), the cross-correlation statistic defined by Equation ( 
written in the Fourier domain as 


161 can be 


Sij = 


f 


Q{f)sKf)sdf)df, 


(Al) 


which for discretely sampled data becomes a weighted sum over 
frequency components, fk. Here the subscripts i and j have 
been introduced to denote particular pulsars — Sij is the cross¬ 
correlation of the data from pulsar i and pulsar j. 

The aim is to choose the filter function Q and a suitable com¬ 
bination of the cross-correlations of different pulsar pairs that is 
‘optimal’ in the sense that it maximises some quantity that is rep¬ 
resentative of the probability of detection. Normally, the optimal 
filter Q is chosen for a given pair of pulsars first and then these 
optimal statistics for each pulsar pair are combined to give the final 
statistic. However, since the dependence of Sij on Q is linear, if 
we consider linear combinations of the S./s these two stages can 
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be done simultaneously and the problem reduces to considering a 
statistic of the form 


X = 2 EE 

k ij 


and finding the optimal combination of the coefficients Xijk- The 
factor of 2 in the preceding equation comes from replacing the in¬ 
tegral over both negative and positive frequencies by a sum over 
positive frequency components only. 

The true optimal statistic will be the combination that max¬ 
imises the DP at a fixed FAP, but as a proxy for this it is usual to 
consider statistics that maximise the S/N, which is the ratio of the 
expected value of a statistic in the presence of a signal, to its 
standard deviation. The standard deviation can either be computed 
in the absence of a signal, (Jq, or in the presence of a signal, cri. 
The A-statistic was constructed to maximise fii/ao. This can be 
thought of as a measure of the significance level at which a given 
source has a good chance of being detected. The noise-only stan¬ 
dard deviation cto determines the FAP, and a given FAP therefore 
corresponds to a certain number of standard deviations, kao- As¬ 
suming a symmetric distribution, a source with /ri/fro = k has a 
50% chance of being detected if the threshold is set to that FAP. 
The B-statistic was constructed to maximise /'^i ■ This is a mea¬ 
sure of how inconsistent such a signal is with noise. The posterior 
distribution in the presence of a signal will have width cri and so 
this statistic measures the number of standard deviations above 0 
the mean of the posterior lies. As we will see below the two statis¬ 
tics are equivalent in the weak-signal regime, but the B-statistic is 
more robust for stronger GWBs. 

The residuals from each pulsar consist of a signal plus addi¬ 
tive noise. The signal is correlated between different pulsar pairs, 
but uncorrelated at different frequencies, while the noise is uncorre¬ 
lated between pulsars and also between frequencies. In the absence 
of a signal we therefore have 

{sijk)o = 0, va.r{sijk)o = Pi{fk)Pjifk) 

COV^Sijk ^ Slmfo ~ 0, COv(Sijfc, SiZTn)o ~ 0 (A3) 

in which Pff) is the power-spectral density of noise in the pulsar 
i. In the presence of a signal we instead have 


{Sijk)\ 

var(sijfc)i 

COv(Sijfc, Slmri)\ 
CGvi^Sijk, Silmf 


t^ijShifk), 

[P{fk)+SH{fk)mifk) + SHifk)] 

prlsUfk) 

[TilTj-m + t'imT jl]Sh{fk)3kn, 
P{fk)rjiSh{fk) 

pf jl p tl'ijT'il]Sk,{fk)3km- (A4) 


Using these results we can compute /ii, (tq and cri for the statistic 
defined by equation (|A2^ 


(A) 

k ij 

(A5) 

var(A)o 

— 2 ^ ^ ^ ^ ^ ^ ^ijk [Eo]ij,im (/fc) Aimfc 

k ij Im 

(A6) 

var(A)i 

— 2 ^ ^ ^ ^ ^ ^ ^ijk (/fc) Aimfc 

(A7) 


k ij Im 


in which T,o{fk) and Ei(/j;) are A^pp x A^pp covariance matrices 
for each fk, where A^pp is the number of pulsar pairs, labelled by 


ij, and 

^Oij,ij{fk') ~ var(sijfc)o, 

^Oij,lm{fk') ~ COV^Sijk, Slmk^O (A8) 

and similarly for Ei. 


A2 A-statistic 

Maximising fj,/ao = (X)/-\/var(A)o is equivalent to maximis¬ 
ing {X} subject to the constraint var(A)o = 1, which straightfor¬ 
wardly yields the result 

Xijk OC ^]ij,lm{fk)t'lmSh{fk)- (A9) 


The matrices Eo(/fc) are diagonal and so we obtain the A-statistic 


^a = 2EE 

k ij 


PjjShoifk) 

Pifk)PAfkf'^'^' 


(AlO) 


where we have introduced Shoifk) to denote the value of the spec¬ 
tral density used to construct the statistic, which is fixed, to distin¬ 
guish it from Sh (fk), the actual value of the spectral density in the 
background, which in general might not be equal to Shoifk)- The 
expected value of the A-statistic in the presence of a signal and its 
variance in the absence and presence of a signal are then straight¬ 
forwardly given by 


„ t'ijShSho 


2 


2EE 


p2 q2 

^ ij‘JhO 

Pp. ’ 
JPiJPj 


(All) 

(A12) 


2 .^^^"jSl,[lP + Sk][PjPSn]Prj,Sl] 

'^1 rppi2 • 

k ij 


where it should be assumed that Sh, Sho and Pi are evaluated at fk, 
but we have suppressed the argument for compactness. If Sho = Sh 
the ‘A’ signal-to-noise ratio, ni/ao, is 


S/Na 


2EE 

k ij 


p2 q2 2 

^ ij'^h 


PP.i 


(A14) 


which is the formula that is most commonly used in the GWB lit¬ 
erature. 


A3 B-statistic 

The equivalent result for the maximisation of /r/(Ti = 
{X) / ^/vsirfXji takes the same form, with Eq replaced by Ep. 
The matrices Ei(/fe) are not diagonal, but the off-diagonal terms 
are quadratic in Sh{fk) and so it is normal to ignore these terms 
and assume they are sub-dominant relative to the diagonal terms. 
This leads to the B-statistic previously given in Equation j20| l, 

Ab = 

2I'ijSho{fk)Sijk 

V ^ WA p Shoifk)]{PAfk) P Shoifk)] P rislffk) ■ 

(A15) 
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The corresponding expressions for the B-statistic are 

2EE' 


Ml 


2 

O-Q = 


TijShSho 


f 

2EE 


lP. + Snom+Sf,o]+r%Sl, 




/ [[P^ + SHom + s^o] + rj,sl,] 


2 _ o V- V- [[P + Snm + Sh] + r^sl 


/ t,- [[J’» + 5M][P.+5,o]+rTS2o]^ 

and for Sho = Sh, the ‘B’ signal-to-noise ratio, /iijai, is 
S/Nb = 

1/2 

1~i2 ^2 

2EE 


■p2 o2 

ij'^h 


f */ 


pp + s4P, + p,] + sl[i + rl] 


(A16) 


(A17) 


(A18) 


(A19) 


This fo rmula is not equivalent to Equation (17) of |Siemens et al.| 
(2013 I because of the term multiplying in the denominator; 
numerically, the difference is negligible. This expression can be 
derived from Section V.A of|Allen & Romano|(|1999|l. 


A4 Comparison 


These two statistics are equivalent at low GWB amplitudes, but can 
deviate as the amplitude increases. If we assume that we always 
adjust the statistic to match the GWB amplitude so that Sho = Sh, 
then (jQ/(Ti tends to 0 as S'), —>■ oo, while /ai tends to a constant 
value. From Equation HI we see that the DP tends to a finite value 
less than 1 as Sh —> oo. For the A-statistic the limiting value of 
Mi/di is 


Mr 

O’! 


2 J2k ^ijlPiPj 


while for the B-statistic it is 
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Both expressions involve sums over all pulsar pairs. Therefore, if 
all the Pi’s are of similar magnitude, the two expressions scale like 
~ Np, and the limiting value of the DP is very close to 
1. However, if one of the P's is much smaller than the others, the 
A-statistic effectively has fewer terms in it and nijoi scales like 
y/Np. The A-statistic is therefore less robust to having an inhomo¬ 
geneous set of pulsars and so we use the B-statistic to characterise 
GWB detection in this analysis. 

We finish with three comments. 


(i) Firstly, the function Sho appearing in these statistics is not 
known a-priori. In practice, this quantity is what we are trying to 
measure with our observations. In the case of PTAs, it is most likely 
that the background will come from a population of merging SB- 
HBs of the type described in this paper and for which the slope 
should be very close to Sh oc and so fixing this as the 

frequency dependence of Sho is likely to be close to optimal. The 
amplitude is not known a priori and does not cancel out of the ex¬ 
pressions for the B-statistic. One reasonable strategy would be to 
compute the amplitude of a marginally detectable background, i.e., 
fix the FAP and compute the background amplitude that yields a 
~ 50% probability of detection, and use this in Sho- Again this 
should be near optimal when a background is first detected and 


will still yield confident detections for louder backgrounds, even 
though it will be sub-optimal at such amplitudes. Other approaches, 
such as using a set of statistics constructed with a range of different 
functions Sho, or an iterative approach in which an estimate of the 
amplitude obtained from the measured data is used to contract a 
new statistic and so on, would also be possible. An investigation of 
these approaches to data analysis is beyond the scope of this paper. 
However, we found that, in practice, the performance of the statis¬ 
tic did not depend strongly on the assumption about Sho and so for 
all the calculations in this paper we assume perfect knowledge of 
the background and set Sho = Sh- 

(ii) Secondly, in deriving the B-statistic we ignored the off- 
diagonal elements of the cross-correlation matrix. The B-statistic 
will be valid in an intermediate signal regime, but will become 
sub-optimal at large background amplitudes due to the omission 
of these terms. This was also discussed in | Allen & Romano| ( |1999| l, 
where it was argued that these terms were always subdominant for 
observations with GBDs. The importance of these terms has not 
been investigated for PTAs. 

(iii) Finally, neither the A-statistic nor the B-statistic is truly 
optimal in the sense of maximising the DP for a fixed FAP. Not 
only have off-diagonal correlations been ignored in deriving the B- 
statistic, but the procedure of maximising hi/oq or fiijoi is not 
equivalent to maximising the DP. Under a Gaussian assumption we 
need to maximise [/ti — fccro]/cri where fc is a constant that depends 
on the chosen FAP. This maximisation can be done, although it is 
more difficult, and we did not consider it here. Simulations indi¬ 
cate that both the A- and B-statistics are nearly optimal in the weak 
signal limit, while the B-statistic continues to be near optimal for 
higher background amplitudes. 
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